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Abstract 

10 

Purpose: To develop an algorithm for real-time volumetric image reconstruction and 3D 
tumor localization based on a single x-ray projection image for lung cancer radiotherapy. 

Methods: Given a set of volumetric images of a patient at N breathing phases as the 
1 5 training data, we perform deformable image registration between a reference phase and 

the other N-l phases, resulting in N-l deformation vector fields (DVFs). These DVFs can 
be represented efficiently by a few eigenvectors and coefficients obtained from principal 
component analysis (PCA). By varying the PCA coefficients, we can generate new 
DVFs, which, when applied on the reference image, lead to new volumetric images. We 
20 then can reconstruct a volumetric image from a single projection image by optimizing the 

PCA coefficients such that its computed projection matches the measured one. The 3D 
location of the tumor can be derived by applying the inverted DVF on its position in the 
reference image. Our algorithm was implemented on graphics processing units (GPUs) to 
achieve real-time efficiency. We generated the training data using a realistic and dynamic 
25 mathematical phantom with 10 breathing phases. The testing data were 360 cone beam 

projections corresponding to one gantry rotation, simulated using the same phantom with 
a 50% increase in breathing amplitude. 

Results: The average relative image intensity error of the reconstructed volumetric 
30 images is 6.9% ± 2.4%. The average 3D tumor localization error is 0.8 mm ± 0.5 mm. On 

an NVIDIA Tesla CI 060 GPU card, the average computation time for reconstructing a 
volumetric image from each projection is 0.24 seconds (range: 0.17 and 0.35 seconds). 

Conclusions: We have shown the feasibility of reconstructing volumetric images and 
35 localizing tumor positions in 3D in near real-time from a single x-ray image. 
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40 1. Introduction 

Management of tumor motion is a challenging and important problem for conformal lung cancer 
radiotherapy. Poorly managed tumor motion can lead to poor target coverage and an 
unnecessarily high dose to normal tissues 1 . Therefore, precise knowledge of real-time lung tumor 

45 motion during the treatment delivery is essential for the effectiveness of lung cancer 
radiotherapy 2 " 8 . Many treatment machines are equipped with an on-board imaging system, which 
usually consists of one kilo-voltage x-ray source and one imager. However, it is very difficult to 
localize lung tumors in 3D space based on a single x-ray projection image, especially without the 
aid of implanted radio-opaque markers. Zeng et al. estimated 3D respiratory motion from cone 

50 beam projections based on a generic B-spline motion model 9 . However, because of the large 
number of parameters in the model, many projections over a 180° rotation have to be used. The 
estimation is a retrospective and lengthy process (computation takes several hours on MATLAB) 
and cannot be used for real-time tumor localization. Therefore, a more efficient lung motion 
model is necessary for such purposes. Zhang et al. developed a lung motion model based on 

55 principal component analysis (PCA), which can efficiently represent the lung motion with only a 
few eigenvectors and coefficients 1 . The PCA lung motion model was recently shown to bear a 
close relationship with the physiological 5D lung motion model proposed by Low et al. 1 1 on a 
theoretical basis 12 . What's more, the implicit regularization of the PCA model allows one to 
derive dynamic lung motion with very limited information, and yet in a reasonably accurate way. 



In this work, we describe a new algorithm based on the PCA lung motion model and demonstrate 
its ability of reconstructing volumetric images and extracting 3D tumor motion information in 
real time from a single x-ray projection in a non-invasive (no marker implantation required), 
accurate, and efficient way. 



2. Methods 

Given a set of volumetric images of a patient at N breathing phases as the training data, we 
perform deformable image registration between a reference phase and the other N-l phases, 

70 resulting in a set of N-l deformation vector fields (DVFs). This set of DVFs can be represented 
efficiently by a few eigenvectors and coefficients obtained from PCA. By varying the PCA 
coefficients, we can generate new DVFs, which, when applied on the reference image, lead to 
new volumetric images. We then can reconstruct a volumetric image from a single projection 
image by optimizing the PCA coefficients such that its computed projection image matches the 

75 measured one. The 3D location of the tumor can be derived by applying the inverted DVF on its 
position in the reference image. In reality, the set of training volumetric images can either be 
4DCT available from treatment simulation or 4DCBCT available during patient setup. 

2.1 PCA lung motion model 



In the PCA model, the DVF relative to a reference image as a function of space and time is 
approximated by a linear combination of the sample mean vector and a few eigenvectors 
corresponding to the largest eigenvalues, i.e., 



85 where u k are the eigenvectors obtained from PCA and are functions of space only. The scalars 

w k {t) are PCA coefficients and are functions of time only. It is worth mentioning that the 

eigenvectors are fixed after PCA and it is the evolution of the PCA coefficients that drives the 
new volumetric image and the dynamic lung motion in real time. 
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90 There are primarily two reasons why the PCA lung motion model is suitable for this work. First, 
PCA provides the best linear representation of the data in the least-mean-square sense 
(efficiency). Second, the PCA motion model imposes inherent regularization on its 
representation. One can show that if two voxels move similarly, their motion represented by PCA 
will also be similar 12 . The combined effect is that a few scalar variables (PCA coefficients) are 

95 sufficient to dynamically deform the lung in a reasonably accurate way. 

2.2 Image reconstruction using PCA model 

After we have obtained a parameterized PCA lung motion model, we seek a set of optimal PCA 
100 coefficients such that the projection of the corresponding volumetric image matches with the 
measured x-ray projection. However, the computed and measured projection images may have 
different intensity levels. Here we assume there exists a linear relationship between them. The 
cost function is: 

min.j(w, a,b) = ||P -f (x,f ) -a- y -b -l|| 2 
s.t. x = x + U • w 

105 where, U is matrix whose columns are the PCA eigenvectors and w is a vector comprised of 
the PCA coefficients to be optimized, x is the parameterized DVF, f is the reference image, f 

is the reconstructed image, y is the measured projection image, and P is a projection matrix 
which computes the projection image of f . For clarity of notation, we have suppressed the time 
index under w , x , y and P . 
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To find the optimal values for w,a,b , the algorithm alternates between the following 2 steps: 

cL/" 

stepl: w, j+1 =w„ -//„•— (3) 

^p2:( a „ +lA+ r =(Y r Y)-Vpf„ +1 (4) 
where, Y = [y,ll, and — = — • — • — = 2-U r • — •P r -(P-f-a-y-6-l) . 



Given the updated PCA coefficients and hence DVF in step 1, the reconstructed image f n+1 is 

found through trilinear interpolation. Accordingly, df/dx has to be consistent with the 

interpolation process in order to get the correct gradient. It turns out that df/dx is a linear 

combination of the spatial gradients of the reference image evaluated at the neighboring eight grid 
120 points, weighted by the appropriate fractional part of the DVF. Step 1 is a gradient descent 
method with variable w and fixed a, b, where the step size fi n is found by Armijo's rule for line 

search. In step 2, the update for a, b is the unique minimizer of the cost function with fixed w . 
Therefore, the cost function always decreases at each step. Note that the cost function is lower 
bounded by zero. The above alternating algorithm is guaranteed to converge for all practical 
125 purposes. The algorithm stops whenever the norm of the gradient is sufficiently small or the 
maximum number of iterations is reached. At this point, we have reconstructed the volumetric 
image which is the best estimate of the current patient geometry in the least mean square sense, 
given the reference image, the parameterized DVF, and the measured projection image. 

130 2.3 Tumor localization via deformation inversion 



In order to get the current tumor position, it is important to distinguish between two different 
kinds of DVFs: push-forward DVF and pull-back DVF. The DVF found by Eq. (2) is a pull-back 
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DVF. It cannot be used directly to calculate the new tumor position. To do that, we need its 
135 inverse, i.e., the push-forward DVF. Here, we adopt an efficient fixed-point algorithm for 
deformation inversion 13 and calculate tumor position. 

2.4 GPU implementation 

140 Recently, general -purpose computing on graphic processing units (GPUs) has offered superb 
efficiency for computationally intensive tasks in medical imaging and therapy applications 14 ' 15 . 
To achieve real-time efficiency, we have implemented our algorithm using compute unified 
device architecture (CUDA) as the programming environment. 

145 3. Materials 

The algorithm was tested using a non-uniform rational B-spline (NURBS) based cardiac-torso 
(NCAT) phantom 16 . This mathematical phantom is based on data from the Visible Human 
Project, and is very flexible, maintaining a high level of anatomical realism (e.g., a beating heart, 

150 detailed bronchial trees, etc.). The respiratory motion was developed based on basic knowledge of 
respiratory mechanics. We generated a dynamic NCAT phantom composed of 10 breathing 
phases as the training data, with a 3D tumor motion magnitude of 1.6 cm and a breathing period 
of 4 seconds. The dimension of each volumetric image is: 256x256x120 (voxel size: 2x2x2.5 
mm 3 ). We used the end of exhale (EOE) phase as the reference image and performed deformable 

155 image registration (DIR) between the EOE phase and all other phases. The DIR algorithm used 
here is a fast demons algorithm implemented on GPU 15 . Then PCA was performed on the nine 
DVFs from DIR and three PCA coefficients and eigenvectors were kept in the PCA lung motion 
model. 

160 For testing purposes, we try to reconstruct a volumetric image and derive the 3D tumor location 
from each of 360 cone beam projections, which are generated from the NCAT phantom and 
uniformly distributed over one full gantry rotation. The phantom has a 4 s breathing period and a 
50% increase in breathing amplitude (3D tumor motion magnitude: 2.4 cm) relative to the 
training data. The gantry rotation lasts one minute, resulting in 15 breathing cycles and 24 

165 projections per cycle. Corresponding to 24 projections, we have 24 volumetric images at 24 
breathing phases. These 24 volumetric images will be used as ground truth test images to evaluate 
the accuracy of the reconstructed images. The imager has a physical size of 40x30 cm 2 . For 
efficiency considerations, we down-sample every measured projection image to a resolution of 
200x150 (pixel size: 2x2 mm 2 ). 

170 

To quantify the accuracy of our reconstruction algorithm, we use the relative image error defined 

as: e = <^ j .[f i -f*^' j^.f* 2 , where / is the reconstructed volumetric image, /* is the ground 

truth image from the NCAT phantom, and the summation is taken over all the voxel indices. We 
quantify the localization accuracy by the 3D root-mean-square (RMS) error. 

175 

The accuracy of our algorithm is limited by the accuracy of training DVFs derived with the 
demons algorithm. The best accuracy we can achieve is the accuracy of the DIR between the 
ground truth test image and the reference image using demons. We therefore deform the reference 
image to the ground truth test image using the same demons algorithm and compute the relative 
180 image error of the deformed image against the ground truth test image. This error, termed as 
deformation error, in contrast to the reconstruction error of the reconstructed image, is used as 
the benchmark to evaluate our algorithm. Similarly, we also compute the 3D RMS deformation 
error. 
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185 4. Results 



Figure 1 shows the "measured" projection image, simulated from the NCAT phantom at the end 
of inhale phase, so that it is maximally different from the reference image. We performed 10 
iterations and further iterations were found to have little influence on the results. The relative 

190 image reconstruction error is initially 35% and approaches to the deformation error (10.9% 
compared with 8.3%) after 9 iterations. Figure 2 shows the sagittal and coronal views of the 
absolute difference images between the reference image and three other images, namely the 
ground truth test image, the image reconstructed by our algorithm, and the deformed reference 
image using demons. The 3D RMS tumor localization error is about 0.9 mm. (anterior-posterior: 

195 0.2 mm; lateral: 0.8 mm; superior-inferior: 0.4 mm). 
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Fig. 1. Left: "measured" projection of the test image at a right posterior oblique (RPO) angle; Right: 
objective value and the relative image reconstruction error at each iteration. 
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Fig. 2. Left column: sagittal view of the absolute difference image between: a) ground truth test and 
reference images, b) ground truth test image and image reconstructed using a single projection, c) ground 
truth test image and the deformed reference image using demons. Tumor is a round object near the center 
of the slice. Right column: same as left column, except for coronal view. Tumor is a round object in the 
right lung. 
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Overall, for the 360 cone beam projections, the average relative image reconstruction error is 
210 6.9% ± 2.4%. The average 3D tumor localization error is 0.8 mm ± 0.5 mm and is not affected by 
projection angles (see Fig. 3). In comparison, the average relative image deformation error is 
5.4% ± 2.2% and 3D RMS tumor localization error from deformation is 0.2 mm± 0.1 mm. 




150 200 
Projection angle 

215 Fig. 3. Top row: relative image error between the ground truth test image and: reference image (red); image 
reconstructed using the proposed algorithm (blue) as a function of cone beam projection angle. Bottom 
row: same as top row, except for 3D localization error. 



We initialized the PCA coefficients to those at its previous frame considering the high-frequency 
220 image acquisition (6 Hz) relative to breathing. The image reconstruction and tumor localization 
for each projection is achieved in less than 0.4 seconds using an NVIDIA Tesla CI 060 GPU card. 
Particularly, depending on the difference between the test image and reference image, the 
algorithm converges between 0.17 seconds and 0.35 seconds (average 0.24 seconds). Compared 
with an implementation on MATLAB 7.7 running on a PC with a Quad 2.67 GHz CPU, which 
225 takes around 15 minutes to converge, our GPU version achieves a speedup factor of more than 
2500. 



5. Discussions and Conclusion 



230 We have shown that it is feasible to reconstruct volumetric images and localize 3D tumor 
positions from a single x-ray projection in a non-invasive, accurate, and efficient way. The 
relative image reconstruction error is only marginally larger than its lower bound defined by the 
deformation error. The average 3D RMS tumor localization error is below 1 mm. By utilizing the 
massive computing power of GPUs, we were able to reconstruct a volumetric image and derive 

235 3D tumor locations from one projection within 0.4 seconds. 

To investigate if the projection image resolution of 200x 150 is sufficient, we tested the algorithm 
on a higher resolution of 1000x750 (pixel size: 0.4x0.4 mm 2 ). The improvement of using the 
fine -resolution projection images is negligible (both image error and localization error are the 
240 same as above). In fact, further reduction in imager resolution does not have noticeable effects on 
the results unless it drops below 40x30, for the phantom case tested here. This demonstrates the 
ability of the PCA model to infer lung motion using limited information. For real clinical cases, 
higher resolution may be needed which will be investigated in our future work. 

245 We plan to test the accuracy of our algorithm on clinical data. There are several complicating 
factors that may affect the accuracy of the current algorithm. First, the x-ray energies used to 
generate the training volumetric images (such as 4DCT) and cone beam projection images may be 
different. Second, there will be scattering and other physical effects that may degrade the quality 
of the cone beam projection images in patient data. In both cases, a linear relationship between 



6 



250 the image intensity of the computed and measured projection images may not be accurate. Some 
preprocessing (e.g., a nonlinear transform) may be needed. We will investigate the application of 
more general similarity measures such as mutual information. This work has focused on x-ray 
projections with rotational geometry. The same principle can be easily applied to those with 
fixed-angle geometry, such as fluoroscopy. The only difference is that the projection matrix will 

255 be constant for fixed-angle geometry instead of varying at each angle for the rotational geometry 
considered in this work. To further speed up the computation, we can also use the predicted PCA 
coefficients from previous histories to serve as the initial guess for optimization. This procedure 
might be important when the imaging frequency is low. 
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